"""
Its easy to create leaflet maps wih Python!
python leaflet_gshhs.py ~/Downloads/gshhg-bin-2.2.2/gshhs_h.b
"""
import sys, math, Image, ImageDraw, os

from gshhs import Parser

LATMIN = 49.9
LATMAX = 58.7
LONMIN = -12
LONMAX = 1.85
ZOOM=12
LAND='lightgray'
SEA='lightgray'
PATH=os.getcwd()+'/tiles'
if not os.path.isdir(PATH):
    os.mkdir(PATH)

im = Image.new('RGB', (256,256), SEA)
draw = ImageDraw.Draw(im)
im.save('%s/SEA.png' % PATH)
im = Image.new('RGB', (256,256), LAND)
draw = ImageDraw.Draw(im)
im.save('%s/LAND.png' % PATH)

class Simul(object):
    """
    Simultaneous equation solver for drawing on arbitary tiles.
    """

    def __init__(self,x,y,a,b):
        """
        mx+c=a
        my+c=b

        m(x-y)=a-b
        c=a-mx
        """
        m=self.m=(a-b)/(x-y)
        c=self.c=a-m*x

    def __call__(self,x):
        """
        """
        return self.m*x+self.c

def deg2num(lat_deg, lon_deg, zoom): # http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

def num2deg(xtile, ytile, zoom): # http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)

class Application(Parser):
    """
    """

    def __init__(self):
        """
        """
        self._tiles = {} # key by x,y,zoom to set of polygone
        self._polygons = {} # key by h.pid
        self._headers = {} # key by h.pid
        Parser.__init__(self, sys.argv[1])

    def header(self, h):
        """
        """
        self._header = None
        if h.c == 'L':
            return
        if h.container>=0:
            return
        self._header = h
        self._polygon = []
        assert h.pid not in self._polygons
        self._polygons[h.pid] = self._polygon
        self._headers[h.pid] = h

    def point(self, p):
        """
        """
        if self._header is None:
            return
        self._polygon.append((p.lat, p.lon))
        if LATMIN<=p.lat<=LATMAX and LONMIN<p.lon<LONMAX:
            for z in range(ZOOM+1):
                if z<6:
                    continue
                x,y=deg2num(p.lat, p.lon, z)
                self._tiles.setdefault((z,x,y), set()).add(self._header.pid)

    def create_tile(self,key):
        """
        """
        z,x,y = key
        zpath = '%s/%s' % (PATH,z)
        if not os.path.isdir(zpath):
            print zpath
            os.mkdir(zpath)
        xpath = '%s/%s' % (zpath, x)
        if not os.path.isdir(xpath):
            print xpath
            os.mkdir(xpath)
        ypath = '%s/%s.png' % (xpath,y)
        if os.path.exists(ypath):
            return
        adjlat1, adjlon1 = num2deg(x,y,z)
        adjlat2, adjlon2 = num2deg(x+1,y+1,z)
        im = Image.new('RGB', (256,256), SEA)
        draw = ImageDraw.Draw(im)
        print 'edge', ypath
        calcx = Simul(adjlon1, adjlon2, 0, 255)
        calcy = Simul(adjlat2, adjlat1, 255, 0)
        for hpid in self._tiles[key]:
            points = []
            for lat,lon in self._polygons[hpid]:
                points.append((calcx(lon), calcy(lat)))
            if SEA==LAND:
                draw.polygon(points, fill=LAND, outline='black')
            else:
                draw.polygon(points, fill=LAND)
        im.save(ypath)

    def fill_tile(self, z,x,y):
        """
        """
        ypath = '%s/%s/%s/%s.png' % (PATH,z,x,y)
        if os.path.exists(ypath):
            return
        if LAND==SEA:
            print 'link', ypath
            os.system('ln -sf %s/LAND.png %s' % (PATH, ypath))
            return
        print 'fill', ypath
        adjlat1, adjlon1 = num2deg(x,y,z)
        adjlat2, adjlon2 = num2deg(x+1,y+1,z)
        im = Image.new('RGB', (256,256), SEA)
        draw = ImageDraw.Draw(im)
        calcx = Simul(adjlon1, adjlon2, 0, 255)
        calcy = Simul(adjlat2, adjlat1, 255, 0)
        for key in sorted(self._tiles):
            _z,_x,_y = key
            if _z!=z:
                continue
            for hpid in self._tiles[key]:
                header = self._headers[hpid]
                if (adjlon1<header.west or adjlon1>header.east) and (adjlat1<header.south or adjlat1>header.north):
                    continue
                points = []
                for lat,lon in self._polygons[hpid]:
                    points.append((calcx(lon), calcy(lat)))
                draw.polygon(points, fill=LAND)
        im.save(ypath)

    def done(self):
        """
        """
        tile_range = {} # key by z to x's and y's
        for key in sorted(self._tiles):
            z,x,y = key
            tile_range.setdefault(z, {})
            tile_range[z].setdefault('x',[]).append(x)
            tile_range[z].setdefault('y',[]).append(y)
            self.create_tile(key)
        print 'tiles', len(self._tiles)
        for z in sorted(tile_range):
            for x in range(min(tile_range[z]['x']), max(tile_range[z]['x'])+1, 1):
                for y in range(min(tile_range[z]['y']), max(tile_range[z]['y'])+1, 1):
                    self.fill_tile(z,x,y)

Application()
